

loess.time.series<-function(data, dv, group, fun, cut){
	library(doBy)
##compute daily hit rate in the data

form<-paste(as.character(dv),"~",as.character(group))
form2<-as.formula(form)

means.all.day<-summaryBy(form2, data=data, FUN=fun, na.rm=TRUE)##get monthly hit rates for weapon related stops
means.all.day<-na.omit(means.all.day)
means.all.day<-means.all.day[order(means.all.day$dayno, decreasing=F), ]


outcome<-names(means.all.day)[2]

##pre-treatment days
dd1<-na.omit(means.all.day[means.all.day$dayno<cut, c(group, outcome)])


##post-treatment days
dd2<-na.omit(means.all.day[means.all.day$dayno>=cut, c(group, outcome )])


set.seed(5103)##set seed to generate reproducible results
##determine spans through 5-fold cross-validation
errors<-crossvalidate.number(dd1, outcome, group, 5 )
best.span1<-errors$spans[errors$means==min(errors$means)]
set.seed(5103)
errors<-crossvalidate.number(dd2, outcome, group, 5 )
best.span2<-errors$spans[errors$means==min(errors$means)]


##estimate loess models and save predictions
form3<-paste(as.character(outcome),"~",as.character(group))
form4<-as.formula(form3)

model1<-predict(loess(form4, data=dd1, span=best.span1, model=TRUE, x=TRUE, y=TRUE), se=T)
fit1<-model1$fit
lb1<-model1$fit-1.96*model1$s
ub1<-model1$fit+1.96*model1$s
i1 <- order(dd1$dayno)



model2<-predict(loess(form4, data=dd2, span=best.span2, model=TRUE, x=TRUE, y=TRUE), se=T)
fit2<-model2$fit
lb2<-model2$fit-1.96*model2$s
ub2<-model2$fit+1.96*model2$s
i2 <- order(dd2$dayno)

lower<-list(fit1, lb1, ub1, i1)
names(lower)<-c("fit1","lb1","ub1","i1")
upper<-list(fit2, lb2, ub2, i2)
names(upper)<-c("fit2","lb2","ub2","i2")
results<-list(lower, upper, means.all.day, dd1, dd2)
names(results)<-c("lower","upper","means.all.day","dd1","dd2")
return(results)
}


